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Using extensive Monte Carlo simulations, we clarify the critical behaviour of the 3 dimensional 
simple cubic Ising Fully Frustrated system. We find two transition temperatures and two long range 
ordered phases. Within the present numerical accuracy, the transition at higher temperature is 
found to be second order and we have extracted the standard critical exponent using finite size 
scaling method. On the other hand, the transition at lower temperature is found to be first order. 
It is argued that entropy plays a major role on determining the low temperature state. 



I. INTRODUCTION 

One of the nowadays trends in statistical physics is the study of disorder. The Spin Glass (SG) 
problem is the most typical subject in the field which, despite huge theoretical and experimental 
works, has not yet been solved. Besides disorder Spin Glass involves also frustration which is known 
to give rise to peculiar, non-trivial, behaviors, especially in Ising systems. It is therefore considered 
as a promising approach to first understand Fully Frustrated (FF) systems without disorder and then 
come closer to Spin Glass by adding disordeJj. From this point of view one of the present authors 
(L.W.B.) and Campbell investigated the critical behavior along the line from the 3 dimensional (3D) 
FF Ising model to the ± J Ising SG model! In order to get further insights along this line, we have 
re-investigatsd,tke 3D FF Ising system on a simple cubic lattice which was once studied more than 
a decade agoau~u but has been left without thorough understandings. The purpose of the present 
paper is to report our new numerical results on nature of the phase transitions in this model. 

Chui et al.u were the first to point out, by means of the Bethe Peierls approach, that the system 
exhibits a first-order phase transition at a finite temperature. They also estimated the degeneracy 

L± 

of ground states which grows proportionally to 2 ■* with L being the linear dimension of the system. 
Blanckschtein et aD (hereafter referred to as BMB) developed a renormalization group (RG) analysis 
(e-expansion) based on the Landau-Ginzburg- Wilson (LGW) Hamiltonian which describes critical 
behavior of the system of interest. Their main result is that the transition is a weakly first order one. 
They also obtained explicitly the 16-fold degenerated ordered configurations of site magnetization 
just below the transition temperature. Grestu attempted to prove BMB's results numerically but 
the result was inconclusive since the signature of a first-order transition, i.e., a bimodal distribution 
in the energy histogram, was not seen. n 

The first numerical study was done by Kirkpatricka which showed that this system exhibit a 
second order phase transition at T c = 1.25J while geometrical argument leads to the absence of 
phase transition. Diep et ata (hereafter referred to as DLN) carried out the most extensive numerical 
study on this system and obtained the following results. The transition from the paramagnetic phase 
to an ordered phase occurs at T c — 1.355J for the infinite size limit and is of second order. At lower 
temperatures the system moves to another phase, in which disorder is located along unidirectional 
lines which in turn form a periodic array. DLN|-|also mentioned the existence of a crossover behavior 
at a temperature higher than T c . Narita et alD studied the system by looking at configurations of 
wrong (or unsatisfied) bonds, and suggested occurrence of a Kosterlitz-Thouless like transition and 
also note the existence of two and possibly three transitions. 

In the present workqby making combined use of the standard Monte Carlo (MC) method and the 
exchange MC methodB, the nature of the critical behavior in the 3D FF Ising system is clarified. 
There exists two ordered phases both with long range ordered magnetization patterns. The pattern 
in the higher temperature phase below the transition temperature T c i is one of the 16-fold degener- 
ated states predicted by the BMB theory, while the one in the lower temperature phase below the 
transition temperature T C 2 is one of the 24-fold degenerated states found by DLN at lowest tem- 
peratures. A mechanism to derive this lower temperature phase, which we call the "4 J-excitation" , 
is proposed. Within the present numerical accuracy, the transition at T c i (T C 2) is of second (first) 
order. The crossover found by DLN above T c \ tends to vanish for larger systems. 
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The paper is organized as follows. In the next section, we will present the model and the methods 
we have used. Geometrical considerations about the ground states and the 4J-excitation are given 
in Sec. Ill . In Sec. ^ the simulated results for various quantities, in the whole temperature range, 



are presented, and in Sec. ^the nature of the two ordered phases and the associated transitions are 
discussed in some details. The last section is devoted to the concluding remarks. 



II. MODEL AND METHODS 



The model studied is described by the standard Hamiltonian, 



H = - 



E 

<i,]> 



J ij SiSj 



where Si's are Ising spins and Jy's are nearest-neighbor interactions which take value +1 or -1 
following patterns described in Fig. |l|.a. We call it the DLN model q Although its bond pattern is 
different from the one of Fig. |l|b (BMB model), the critical behaviour of the two models is expected 
to be the same because every plaquette is frustrated in both systems and, in fact, one can go from one 
system to the other by gauge transformation. However in the present simulation we have examined 
only the DLN model. 





FIG. 1. Structure of the FF 3D SC (a) DLN and (b) Comb-type (BMB) model. The (wavy) thin line represent the 
(anti-)ferromagnetic bond. The numbers labels the sub-lattices used during the simulation. 

In order to investigate critical properties of the model, we have used two Monte Carlo (MC) 
techniques both making use of the heat-bath updating. One is the standard MC method which 
has been used to obtain time evolution, including gradual cooling or heating process, of quantities 
of interest. The other is the exchange MC methods. The basic idea of this method is to run 
simultaneously several replicas having a common bond configuration but attached to heat baths 
of different temperatures, and to exchange stochastically a pair of them according to the detailed 
balance condition for the combined system of all the replicas. The exchange MC (EMC) method is 
quite efficient to examine equilibrium properties of such a system with many degenerate states. 

In the present work the following quantities are evaluated for systems with N — L 3 spins under 
the periodic boundary condition: 

(o) Total internal energy, E, and specific heat C v . 

(b) Overlap distribution defined by 



iv mcB / JV 



(2) 



where 7V mcs is the number of MC steps, and a, b label two different replicas. 



1 Although the first study was done by Kirkpatrick, we will rely on the paper by Diep et al. and hence the name DLN 
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(c) Sub- lattice energy, e a , defined by 

iV mcB n.n. 

e - = |f E iv~: E E -JijSiimit), (3) 

i <E t — 1 j 

where a indicates the sub-lattices (a = 1,...,8) which are defined by the eight corners of one of the unit cube (see 
Fig. ELa). 



(d) Site magnetization defined by 



M * = N~ E *(*)• (4) 



III. GROUND STATES AND 4J-EXCITATION 

A ground state of the Ising system described by eq. |l| is a spin configuration having the minimum 
number of wrong (or unsatisfied) bonds which are defined as —JijSiSj = +J. Since all plaquettes 
in the present 3D FF model are frustrated, a ground state is obtained when all the plaquettes have 
the minimum number of wrong bonds which is unity. For a ground state of the unit cube, it is easy 
to see that only three wrong bonds are required, thereby none of them are on a common plaquette 
and there are two corners which are not touched by them. Actually there are eight (multiplied by 
two when global inversion of spin is taken into account) ways of minimizing the energy of the unit 
cube as shown in Fig. The generic way to create a ground state of the total system, is to 




e f g h 

FIG. 2. The eight ways of minimizing the energy of the unit cube with three wrong bonds (represented by the bold lines). 

pile up these cubes with the only constraint that the wrong bond is at the same place for each 
common plaquette. For example, the neighbor of a in the x direction can be either c or g. Let us 
introduce what we will call periodic ground states. If one take the c cube as a neighbor of a, the 
wrong bond oriented in the x direction is continuous (in the g case it is not) and the next neighbor 
of c in the x direction can be a. So the pattern can be acac... In this case the position of the 
wrong bond perpendicular to the common plaquette has been taking into account. There are only 
sixteen periodic ground states since once the first cube is set there is no choice for the neighboring 
cube (which gives a degeneracy of 8 times 2 for the spin reversal symmetry). In these periodic 
ground states only two kinds of sites exist when one looks at the sub-lattice energy e a ; 6/8 sites 
with e a = — 2J and 2/8 sites with ep — — 6J. 
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FIG. 3. ALL the spins of the line AA, BB or CC can be inverted without changing the total energy 

As first pointed out by KirkpatrickS, in each of the periodic ground states, there are lines of 
spins that can be inverted without changing the total energy of the system. They are the lines AA, 
BB, and CC in Fig. ||.a and were called by linear-chain excitations by DLN. Two wrong bonds are 
connected to each sites of these lines. The choice of the lines inverted is restricted to be in the 
same direction. For example if one inverts the AA line, one site on the adjacent CC line has only 
one wrong bond and thus inversion of this spin (or this CC line) costs excess energy. Therefore one 
needs to count how much lines can be inverted once the direction is chosen. That number is simply 

L 2 

given by L 2 /A and so degeneracy of the ground states is 2~ , which is the same order of magnitude 
as Chui et afl have found. One of the consequence of this inversion is that the neighboring sites of 
this line which had energy —23 and — 6J before the inversion now have — 4J. 

A further interesting fact is that there is an actual process to invert all spins on one line which 
costs an energy of only 4J. In fact if one spin is flipped, for example, the site i on the BB line in 
Fig. ^.a, the wrong-bond pattern changes to the one shown in Fig. ^b, which has higher energy by 
4J before the flip. But now internal field at the two neighboring sites vanishes so that spins on the 
sites can flip without cost of energy. This mechanism was called kink-anti-kink pair by Kirkpatricku 
but here we call it a 4 J-excitation. With the aid of one 4 J-excitation all spins on the BB line can be 
inverted. Strictly speaking, this holds under the periodic boundary condition we have adopted. It is 
important to notice that the presence of one 4J-excitation is associated with entropy of L In 2 which 
is much larger than its energy 4J. This indicates, and will be demonstrated below that entropy 
plays an important role on determining an ordered phase, if any, at lowest temperatures. 



IV. RESULTS 



In this section we present the data in the whole temperature range simulated. They strongly 
suggest the existence of two phase transitions in the present FF model. The details of their critical 
nature will be discussed in the next section. 



A. Specific heat and energy 

The results of the specific heat and the energy simulated by the exchange MC method are shown 
in Fig. ^. They indicate a second order phase transition at the temperature T c i(~ 1.35 J) and 
existence of another transition at T C 2(— 0.70J). 
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FIG. 4. The (•) represent the results of the fluctuation of the energy while the solid lines represents the derivative of the 
energy. The (o) represent the energy. L=24 



B. Overlap distribution 

The overlap distribution P(q) simulated by the exchange MC method is shown in Fig. [|. It 
exhibits quite different shapes either from the one of a simple ferromagnet (with only double peaks 
at q = ±m 2 , m being the uniform magnetization) or from the one of mean- field SG model (with 
double peaks at q = ±f/EA, 9ea being the Edwards- Anderson order parameter, and continuous 
weight between the peaks). The characteristic features of P(q) of the present FF model are as 
follows: 

1. It has a single peak centered at q = in the paramagnetic phase (T > T c i). 

2. At T c i > T > T C 2, five peaks appear. The peaks are narrowed when the size of the system increases. This indicates that 
there exist several states which are thermodynamically stable. We call the phase of this temperature range the high-T 
phase. 

3. At T = 0.75 ~ T C 2, P(q) exhibits a peak at q = with a continuous distribution in both sides of the peak. This implies 
a certain disordered configuration. 

4. Below T C 2 a seven-peak structure appears, indicating occurrence of another ordered state. We call the phase of this 
temperature range the low-T phase. 
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FIG. 5. Overlap Distribution for the FF system. L=24 



C. Sub-lattice and site quantities 

We have performed standard MC simulations of gradual cooling or heating processes and have 
looked at evolution of the sub-lattice energy e a . The result are shown in Fig. f| The lines represent 
evolution with temperature of the eight e a . For each point 10 4 MCS of annealing were used and the 
next 10 4 MCS were used to get e Q . Each curves in the figure represent the result of a single MC 
run but we have checked that other runs give qualitatively the same results. 




FIG. 6. Typical sub-lattice Energy, (a) is obtained by slow cooling while (b) is obtain during the heating of a ferromagnetic 
spin configuration from T = 0.50. All the data points were obtained during the same run for a system of size L = 24 

From Fig. ^.a we can extract the following characteristic features of the gradually cooling process. 
First of all the sub-lattices are linked two by two. In the high-T phase (T c i > T > T C 2) there exist 
two kinds of sub-lattices. One with lower energy, ei^'s, consists of two sub-lattices which are 



directly checked to be located at opposite corners of the cube. 



e« 's, consists of the other six sub-lattices. As the system is cooled gradually ei L ' 



The other with higher energies, 
s decrease, while 



G 



e„ s stay nearly constant (~ —2 J). In this phase the system thus exhibits a symmetry similar to 
the one of the periodic ground state mentioned previously. In the low-T phase (T C 2 > T) there seem 
to exist three kinds of sub-lattices. Four sub-lattices have e a of — 2J. The four other sub-lattices 
are grouped two by two with two different e a which are symmetric around — 4J. This symmetry 
has already been discussed in the discussion of the periodic ground state when the inversion of lines 
of spin is introduced. But this symmetry doesn't hold for case b, where the system was heated 
from the ferromagnetic ground state, because the annealing time was not long enough at this low 
temperature. But it implies that the ferromagnetic configuration is not stable at finite temperatures. 
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FIG. 7. Histogram of the site-magnetization 



By the same cooling procedure as above described, we have examined site magnetization, {Mf }. 
Its distributions in magnitude at four typical temperatures are shown in Fig. Q At T — 1.10, in the 
high-T phase, there are four peaks. However we discuss here about distributions against absolute 
magnitude of {Mf }, since their signs depend on configurations simulated, similarly to those of the 
sixteen periodic ground states. In this sense there are two kinds of {Mf}, one with smaller \Mi\ 
and the other with larger one. The integration of the peaks centered around \mi\ and | rri2 1 yields 
6/8 and 2/8, respectively. Actually, when the corresponding spatial pattern of {Mf} is visualized, 
it is almost perfectly periodic, and agrees with the one derived by the BMB (see Section 5.1). When 
temperature approaches to T C 2 the position of the two peaks tend to saturate to the limiting value, 
i.e., unity. 

At T = 0.70J and 0.60J, in the low-T phase, there is a flat bump around zero. The integration 
of this bump yields weight of 2/8, while the two sharp peaks at ±1 have weight of 6/8. In the 
corresponding spatial patterns of {Mf} we have observed that lines with nearly zero magnetization 
form a periodic lattice, as found by DLN. This configuration is attributed to the 4J-excitation 
mentioned previousely (see Section 5.3). 



V. DISCUSSIONS 



The results of our MC simulations presented in the previous section reveal that there exist two 
ordered phases in the DLN model. In the following we discuss the nature of the two phase transitions 
in detail. 
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A. High-T phase 



When the system is cooled from the paramagnetic phase, an ordered phase appears at T = T c i — 
1.35 J (see Fig. §), which we have called the high-T phase. It is in good agreement with the one 
theoretically predicted by BMB, which can be easily gauge-transformed from their comb- type model 
to the one used by DLN. 

In the DLN model the magnetic unit cell consists of eight unit cubes shown in Fig. ^. In the BMB 
model, on the other hand, it consists of sixteen cubes. There exist sixteen degenerated equilibrium 
configurations including the spin up-down symmetry. Each elementary cube has \M!*\ ~ m at two 
diagonally opposite sites, and \M*\ ~ m/a with a > 1 at other sites with signs specified in such a 
way that three wrong bonds do not intersect and that they do not touch the two sites with larger 
\M"\. We note that these configurations are different from the periodic grouad states for which 
= 1, and that they minimize the LGW Hamiltonian of the present FF systemu. 

The ordered configurations mentioned above are quite consistent with our simulated data de- 
scribed in the previous section. Concerned with the overlap distribution P(q) at T = 1.07 shown 
in Fig. H the positions and weights of the five peaks are what are expected if the system visits 
the sixteen degenerated configurations with equal probability in the exchange MC simulation. The 
occurrence of these ordered states is confirmed also by the {M* }-histogram at T = 1.10 shown in 
Fig. ft, and by the direct inspection on the spatial pattern of {Mf} (not shown). Furthermore in 
Fig. p we plot the ratio of positions of the two peaks in Fig. fj] against temperature. The results are 
consistent with the theoretical prediction a = \/3o expected to hold at temperature close to T c i- 
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FIG. 8. Ratio of the larger to the smaller site-magnetizations. 

In contrast to the previous arguments by Narita et all and DLN, our simulated data strongly 
suggest that each of the sixteen degenerated configurations is long-range ordered (or thermodynam- 
ically stable). To confirm this we have computed the sub-lattice switching time, r sc . By inspection 
of time evolution of each sub-lattice energy e a at a fixed temperature close to T c i we have observed 
that it fluctuates between the higher and lower branches of e a in Fig. ^ with an average interval 
of r sc which grows with the size of the system. To extract r sc more accurately we have used the 
autocorrelation function q(t) which almost saturates to an equilibrium value before exhibiting an 
exponential decay. This exponential decay, at the later stage, is attributed to the sub-lattice switch- 
ing and its relaxation time is examined at T = 1.32 < T c i) for sizes L varying from 12 to 28. The 
result shown in Fig. indicates that r sc is exponential with an argument growing roughly linearly 
with L. This implies that in this 3D FF model the domain-wall free energy is roughly proportional 
to L (and not to L 2 as for ordinary 3D ferromagnet), and that the sub- lattice switching does not 
occur in the thermodynamic limit. 
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FIG. 9. Time of sub-lattice switching vs the size of the system 

|The technique usually used to know the domain-wall free-energy, is the defect-free-energy analy- 
sis! The principle is to introduce a domain-wall in a system and to see how the energy is affected 
by this domain. If the energy increases with the size than the system exhibit a phase transition. To 
obtain the domain wall free energy one compute 

AF = F A p - F P 

where Fap and Fp axe the free energy of the system for the anti-periodic (AP) boundaries and 
periodic (P) boundaries conditions, respectively. With standard MC simulations the free energy can 
be obtained by : 

/3F(T) = / dp'E{l3'). 

JPmi,n 

In the simulation the fimin is set to 1/10 J and we assume that AF is negligible at f3 m in- The 
results are shown in Figure [l(] where we can see that the maximum of AF also have a nearly linear 
dependence in L. We can also note that the two low temperature phase have a different slopes for 
AF against T. 




FIG. 10. Results for the defect-free-energy analysis for sizes L = 8, 12, 16and20 



B. Critical nature of the transition at T& 




FIG. 11. Size-dependence of the specific heat near T c \. 



According to BMB, the LGW Hamiltonian for describing the phase transition in the present FF 
model is an n = 4 "Heisenberg" model with symmetry breaking terms arising from frustration. 
The RG analysis based on the e-expansion method predicts a weakly first order phase transition. 
Within our present numerical accuracy, however, we have not detected any evidence which supports 
a first-order transition at T c \ (see for example Fig. ^ and further discussions in Section 6). The 
present results are in agreement with the previous numerical works by DLN and Gresttl. 

Assuming a second order transition at T c i = 1.347 ± 0.001, we have obtained the following critical 
exponents: 

a = 0.32 ±0.02, v = 0.56 ± 0.02, /3 = 0.25 ± 0.02, 77 = -0.1 ± 0.02. 

These exponents are different from those of the 0(4) one and we don't know to which class of 
universality they belongs. The exponent a of specific heat is estimated either directly as shown in 
or by the finite-size scaling shown in Fig. |l2| In relation to our results demonstrated in 
it is noted that the apparent crossover within the paramagnetic phase observed by DLN is 
considered to be an artifact due to finite-size effect. Correspondingly, our exponents a and v are 
compatible to those obtained by DLN at temperatures above the crossover whose data are considered 
not affected by the finite-size effect. The r\ value given by DLN (~ 0.28) being obtained at T c i is 
considered to be affected by finite-size effect. 
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FIG. 12. Scaling plot of the Specific heat. 
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The Binder cumulantB method is frequently used to determine the exponents v and /3. From 
the later the exponent r\ is calculated by the scaling relation — = d — 2 + 77. However, for the 
present FF model with the degenerated ordered states below T c i and with the presence of the peak 
at q — in P(q), a direct application of the method is rather complicated. Tfais problem for the 
Binder cumulant had already been pointed out for th« three states potts modelstU. Instead we have 
analyzed the finite size scaling of the following ratictJ. 

2 

i? 2 (L l5 L 2 ) log(l%2±)/Iog(yi) = -2^ +g(L 1 2 / »(T- T c )) 

with g(x) being a scaling function and L2 > L\ 2> L2 — L\ . The factor 2 in the r.h.s. of the above 
equation come out because < q 2 > is given by the following scaling form: 

/ 7\r \ 2 at 



N 



where < .. >t is the thermal average and f(x) another scaling function. Here Gij is the correlation 
function < SiSj >t which is given by 



Gij — G(r = \xi — Xj\) oc r 



-(d-2 + rj) 



exp 



Thus the crossing point of i?2(Li, L2) with different set of {L\,L-2) will be (T c , — ^). The scaling 
plot of Rq(Li, L2) giving Td, v and /3 is shown in Fig. 



-0.5 

-1 
-1.5 

-2 
-2.5 

-3 
-3.5 

-4 



lU.Li] 


[8.121 


1 — t? — 1 


12,16] 




[16,20] 


--- -<>----' 


20,24] 


♦ 



T c = 1.3476 
|1/V=J.724 
« 2 *=-1.72 





(T-T C )L V 



FIG. 13. Scaling plot of R 2 {L 1 ,L 2 ). 



C. Low-T phase and nature of the transition at T — T C 2 



As temperature decreases in the high-T phase, the larger \M" \ become almost saturated to unity. 
Then the arguments based on the LGW Hamiltonian, in which no restriction on magnitudes of Mf 
is imposed, are expected to break down. In fact, as already described in Section [^j, the transition 
to the low-T phase occurs around T = T C 2 — 0.7. 

In each of spin configurations in the low-T phase realized by gradual cooling MC process, one 
fourth of all the chains in one of the three directions are disordered (\M" | — 0), while the others have 
almost saturated value (\M£ | ~ 1). Furthermore these disordered chains form a periodic 2D array. 
Thus the low-T phase is 24- fold degenerated as pointed out by DLN. This degeneracy is responsible 
for the positions and weights of the seven peaks structure in P(q) shown in Fig. ^ We argue that 
the above mentioned order in the low-T phase is attributed to the 4J-excitation introduced at the 
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end of Section fl\ As pointed out there, the presence of one 4J-excitation lowers free energy by 
— TLln2 + 4J as compared with the state without it. This holds true so long as the eight chains 
surrounding the disordered chain (with the 4 J- excitation) are firm and play a role of a cage of the 
latter. Also we can at least check that two 4 J-excitations on two chains perpendicular to each other 
cannot cross freely. Therefore for the system to have maximum free-energy gain a periodic array of 
L /4 disordered chains is realized. 

It is rather hard to check the thermodynamic stability of each of 24 degenerated states in this 
phase, not only because MC dynamics becomes slower at these low temperatures, but also because 
system sizes have to be large enough for an almost ideal random walk of 4 J-excitations to be realized. 
Also free energy barrier between the different states become small as temperature decreases. In fact, 
as seen in Fig. [Uj, magnitudes of defect free energy tend to decrease from around T C 2 with decreasing 
temperature. But since those with the larger L have still the larger values, we expect that each 24 
states are thermodynamically stable. In this context, it is noted that defect free energy at T — 0, is 
easily checked to vanish. This, combined with existence of the 4 J- excitation, implies that stability 
of (or free-energy barriers between) the ordered states so far discussed, is guaranteed solely by the 
entropy effect. 

Lastly we discuss about the nature of the transition at T C 2 between the high-T and low-T phases. 
The transition is expected to be of first order, since it occurs between spin configurations having 
different symmetries with 16-fold and 24-fold degeneracy. Among our numerical data, a peculiar 
shape of P(q) at T = 0.75 shown in Fig. |^, the appearance of a significantly large hysteresis in 
various quantities, and, most plausibly, the energy histogram with double peak structure shown in 
Fig. Fil support this argument. 
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T 

FIG. 14. Histogram of the energy for L — 48 at T = 0.716. The histogram show a double peak structure implying a first 
order phase transition 



VI. CONCLUDING REMARKS 

We have shown by extensive MC simulations nature of the phase transitions in the 3D FF Ising 
model. The driving mechanism of the lower transition at T C 2 agrees with what DLN proposed, 
though we have further introduced the "4J-excitation" mechanism explicitly. We have also claimed 
that this transition is of first order. The nature of the high-T phase below T c i is quite consistent 
with the prediction by the BMB theory. Although our results simulated in the present work strongly 
suggest a second-order phase transition at T c i , we cannot exclude a possibility that it is of first order 
when studied in systems with larger sizes and ati closer temperatures to T c as is the case for other 
systems exhibiting a weak first-order transitioned. The search of this possibility is now underway. 

The present work reveals that the transition behavior of the present FF model are dominated by 
a subtle balance of the entropy effects. They are expected to be affected sensitively if disorder is 
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added. In particular, behavior of the 4J-excitation in presence of small disorder is of interest. The 
analysis along this direction is also underway. 
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